Properties of Vacancy Formation in 4 He hep Crystals at Zero Temperature 

and Fixed Pressure 
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Equation of state of 4 He hep crystals with vacancies is determined at zero temperature using the 
diffusion Monte Carlo technique, an exact ground state zero-temperature method. This allows us 
to extract the formation enthalpy and isobaric formation energy of a single vacancy in otherwise 
perfect helium solid. Results were obtained for pressures up to 160 bar. The isobaric formation 
energy is found to reach a minimum near 57 bar where it is equal to 10.5 ± 1.2 K. At the same 
pressure, the vacancy formation volume exhibits a maximum and reaches the volume of the unit 
cell. This pressure coincides with the pressure interval over which a peak in the supersolid fraction 
of 4 He was observed in a recent experiment. 
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Properties of solid 4 He have regained interest since the 
discovery that it exhibits nonclassical moment of inertia 
[l| and several accompanying phenomena (for a review, 
see Refs 0t3|)- Recently, Kim and Chan studied samples 
grown at fixed pressure [5|. They found that the super- 
solid fraction changed threefold over the studied pressure 
range (up to 135 bar) and had a distinct maximum at 
pressures near 55 bar, which is well above melting pres- 
sure. These results call for a study of possible superso- 
lidity mechanisms under isobaric conditions. 

Crystal defects are believed to be indispensable in the 
mechanisms behind supersolidity [2|, |6(. In particular, 
special attention has been paid in the past to vacancy de- 
fects. Unbound vacancies, whether intrinsic to solid he- 
lium or introduced by experimental conditions, are likely 
to cause supersolidity [J, |7[ . Experimental results for en- 
ergy of vacancy formation in pure 4 He are available from 
the work of Fraass et al. |8J|, with additional data and 
analysis by Blackburn and colleagues [9[ . Numerous the- 
oretical works describe vacancies in solid helium [3, llOr - 
1 181 ] . Results for vacancy formation energy are available 
from finite-temperature calculations using path integral 
Monte Carlo methods LL4|, |l5f and from ground-state cal- 



culations with variational [1 1M13L 1 17 1 , shadow path inte- 
gral ground state [17| and diffusion [7j Monte Carlo meth- 
ods. These calculations agree that individual vacancies 
cost too much energy to exist in the ground state of solid 
helium. However, properties of isobaric vacancy forma- 
tion in 4 He have not yet been reported. Additionally, 
several important vacancy properties, such as vacancy 
volume, are only sensible if considered at fixed pressure. 
To address these concerns, we performed a calculation 
of the equation of state of solid 4 He with the number of 
atoms incommensurate with the number of sites of an hep 
lattice filling the volume. Incommensurate in this sense, 
the solid remained crystalline and had one vacancy de- 



fect per (periodic) simulation volume. Obtaining such an 
equation of state allowed us to extract the dependence of 
vacancy properties computed at both fixed pressure and 
fixed density conditions. Importantly, calculated forma- 
tion energy at fixed pressure is considerably smaller than 
the formation energy at fixed density. It also turned out 
to exhibit a minimum around the same pressure where a 
peak in the supersolid response was observed experimen- 
tally. 

The most important thermodynamic quantity charac- 
terizing the presence of vacancies in fixed-pressure sys- 
tems is the Gibbs free energy [19j, which at zero tem- 
perature reduces to enthalpy Consider the enthalpy 
H(N, 0, P) of a commensurate system consisting of N 
particles at pressure P, and H(N — 1, 1, P), the enthalpy 
of an incommensurate system with TV — 1 particles and 
one vacancy at the same pressure P. If h is the en- 
thalpy per particle and h vac the vacancy formation en- 
thalpy, then H(N, 0, P) = Nh(P) and H(N - 1, 1, P) = 
(N — l)h(P) + /ivac(-P). The vacancy formation enthalpy 
in a solid can therefore be expressed [20( as 

TV- 1 
^vac(P) = H(N-1,1,P)-——H(N,0,P) 
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where p mc and p com are the densities of the incommensu- 
rate and commensurate systems at pressure P. Pressure 
can obtained from the derivative of the energy with re- 
spect to the density p as 

2 0E/N 



P = P Z 



dp 



(2) 



thus allowing us to obtain E(P) and p(P) from the E(p) 
dependence available from the calculations. 



The isobaric vacancy formation energy AEp is the en- 
ergy cost of moving a single atom away from a lattice 
site, while keeping the pressure fixed, 

N - 1 
AE P = E(N- 1,1, P)-^ r E(N 7 0,P). (3) 

AEp determines the formation enthalpy, with the addi- 
tion of the work necessary to free the vacancy volume 
^vac- v i a ^vac = AEp + PQ vac . The incommensurate 
system has an excess of volume with regard to its num- 
ber of atoms. It can be written as V(N — 1,1, P) = 
(N — l)v + f2 V ac where v is the volume normally occu- 
pied per particle, V(N,0,P) — Nv. Hence the vacancy 
volume can be obtained from 

TV- 1 

n vac = V(N-l,l 7 P)-^ r V(N,0,P). (4) 

In previously published numerical simulations, the for- 
mation energy was calculated at fixed density as 

N- 1 



AE p = E{N-l,l,p) 



N 



-E(N,Q,p), 



(5) 



which compares two systems with volume adjusted to 
provide equal densities. While AE p matches AEp at zero 
pressure [2l| . they generally differ at non-zero pressures 
necessary to solidify helium. In fact, AE P slowly ap- 
proaches /i V ac in the thermodynamic limit: h vac — AE P oc 
1/N |22[. Attempting to obtain h vac from AE P is thus 
subject to finite-size effects even if h vac itself converges 
quickly. 

Calculations were made with diffusion Monte Carlo 
(dmc), a statistically exact ground-state method. DMC 
has proven indispensable in understanding properties of 
supcrfluid 4 Hc, andyields excellent results for the prop- 
erties of solid 4 He [7|, [23j , in particular for the equation 
of state. DMC projects the excited states from the ini- 
tial conditions (/> : by advancing through imaginary time 
r a function / = <p G exp[— (H — Eq)t]<Pi, where H is the 
Hamiltonian, Eq is a reference energy, and </> G is an impor- 
tance sampling function, also called gui ding wavefunc- 
tion. For details of the method, see Refs [24l426| . Inter- 
actions were modelled with the Aziz hfd-b(he) potential 
[27l | . While this potential is known to produce a slight 
systematic bias in energy of under 70 mK, it allows for a 
good reproduction of the equation of state P(p) [23l I25IJ . 

A good-quality guiding wavefunction <fi G is necessary to 
efficiently sample the energy from the ground state itself. 
Such a wavefunction has to take into account two-body 
correlations between atoms, to provide lattice symme- 
try, to allow for hopping between lattice sites, and finally 
it needs to satisfy the Bose symmetry. We used a sym- 
metrized form of the Nosanow-Jastrow [28[ wavefunction 
that was recently developed for quantum solids [7J and 
has the form 
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FIG. 1: (Color online) Energy per particle of the commensu- 
rate system (()) and the incommensurate system with one of 
the 180 sites being vacant (o). Error bars are smaller than the 
symbols and are shown within them. Solid lines in the main 
plot area are third-degree polynomial fits as described in the 
text. Lower inset separately shows statistical errors with solid 
symbols. Top inset shows pressures for both systems (pressure 
of the incommensurate solid is higher). Dashed line marks 
melting pressure, below which both systems are metastable 
but remain crystalline. 



where N p and N s are respectively the number of atoms 
and lattice sites, f(r) is the pair correlation function 
and g(r) is a function that localizes atoms to the lat- 
tice sites. \t{ — lk\ denotes distance from an atom 
with index i to the lattice site k. SNJ is both Bose- 
symmetric and provides an excellent spatial order [7j. 
We used g(r) = exp[— l/27r 2 ] and the McMillan form 
f(r) = cxp[— 1/2 (b/r) 5 ]. The parameters were obtained 
from variational optimization of the Nosanow-Jastrow 
wavefunction, resulting in b — 1.12a, where a = 2.556A, 
and 7 = (-6. 120+27. 69 per 3 ) / ^/a . Stability was observed 
for changes of timestep, population size and wavefunction 
parameters over the entire range of densities. 

The commensurate solid was simulated by considering 
a periodic system with 180 atoms in a volume containing 
180 hep lattice sites. This number allows for a nearly 
cubic geometry, with ratio of sides of the simulation box 
equal to 1.06 : 1.02 : 1. The incommensurate crystal was 
simulated by using one less atom while keeping the lat- 
tice intact. Size effects resulting from finite simulation 
box size were compensated by separate variational cal- 
culations involving up to 1440 atoms with the methodol- 
ogy of Ref. [29( but with the inclusion of (1/iV) 2 terms. 
Vacancy images created by periodic boundary conditions 
arc sufficiently separated to neglect their interaction [3fj| . 

The results for the equation of state of both commensu- 
rate and incommensurate systems are presented in Fig. [T] 
The resulting density dependence of energy per particle 
e(p) could be fitted accurately with third-degree polyno- 
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FIG. 2: (Color online) Formation energy of a vacancy at fixed 
density, as defined by Eq. (JSJ. Energy is shown with the solid 
line while the area within one standard deviation is shaded. 
The straight line shows the fit at high densities. Symbols 
show numerical data from Refs []jj| (o), [TJ] (0), Qj] (v), 
[13 (□), El (.) and (♦). 



FIG. 3: (Color online) Vacancy formation enthalpy (red, up- 
per curve) and isobaric formation energy (green, lower curve). 
Shaded areas mark one standard deviation. Dashed line 
shows melting pressure. Diamonds show experimental mea- 
surements from Ref. [3| (see text). 



mials, which in this case are equivalent to the form 

e = e Q + b (p/p - l) 2 + c (p/p - l) 3 (7) 

with the coefficients as follows: for the commensurate 
system, po = 25.15 nm~ 3 , eo = —6.45 K, b = 17.3 K, c = 
15.7 K; for the system with a vacancy, po = 25.51 nm -3 , 
e = -6.34 K, b = 22.5 K, c = 11.7 K. From the depen- 
dence of energy on density, we can extract the formation 
properties according to Eqs fUE]). To accurately calcu- 
late the standard deviation in the quantities calculated 
from e(p), we studied the distribution of each such quan- 
tity as determined directly by the statistical errors in the 
calculation of the energy values. 

Vacancy formation energy at fixed density, computed 
according to Eq. ([5]), is shown in Fig. [2] At large densi- 
ties, AE p grows nearly linearly The apparent bending at 
highest considered densities is not discernible within the 
current error levels. For densities p above 30 nm~ 3 , AE p 
can be approximated within 0.5 K as AE p = k(p/p r — 1) 
with k = 108.2 K and p r = 25.16 nm -3 . It is notable 
that AEp approaches a plateau at low densities. The 
lowest value AE p = 16.5 ± 1.4 K is reached at the lowest 
studied density of p = 27.6 nm -3 . Fixed-density forma- 
tion energies obtained by other groups are also plotted 
for comparison in Fig. Our values are higher over a 
range of densities close to melting. 

By matching the pressures, we were able to compute 
the vacancy formation enthalpy and the isobaric forma- 
tion energy, shown in Fig. [3] The experimental results of 
Fraass et al. [8j are shown in the same figure. Samples in 
that experiment were grown at constant pressure, but the 
experiment itself was reported to be performed at fixed 
volume. The experimental results in Fig. [3] are shown 



at the pressure at which Fraass et al. reports the end of 
the solidification process. The coincidence between the 
isobaric formation energy and the experimental results 
compared in this way calls for more careful examination 
of the thermodynamic conditions under which such ex- 
periments are carried out. 

The minimum that we find in AEp is unexpected and 
is not observed in classical solids [20|. Expanding the 
equation of state of the system with one vacancy around 
density p, one obtains 



As app rox(jp) = AEp(p) _ {N _ 1)/3iBc p*P +(D f}_ 



where AP = Pi nc (p) — P C om(p) is the difference in pres- 
sure between the two systems, /3i nc is the compressibility 
of the incommensurate system and P = P com {p)- Us- 
ing the quantities obtained from our calculations, we find 
that AE P pprox matches the result for AEp to within 1 K, 
including the minimum structure. The minimum in AEp 
occurs around the pressure of 57 bar, where the energy 
lowers to 10.5 ± 1.2 K. This pressure is coincidentally 
rather close to the value where Kim and Chan [5j ob- 
served a maximum in the supersolid fraction. Unlike for 
energy, there is not an extremum in vacancy formation 
enthalpy. The minimum in energy is balanced by the 
character of the dependence of vacancy volume on pres- 
sure. Nonetheless, the second derivative of h vac changes 
sign around this pressure. 

The vacancy formation volume f2 vac , shown in Fig. |4l 
is also non-monotonic and reaches a maximum at 57 bar, 
while its volume in relation to the unit cell volume of 
the incommensurate system V{ a t peaks at 63 bar (see 
left inset in Fig. S) statistical uncertainty does not al- 
low us to distinguish whether f2 vac indeed peaks above 
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FIG. 4: (Color online) Vacancy formation volume as a func- 
tion of pressure [Eq. Q]. Shaded area indicates statistical 
error. The curve peaks at P = 57 bar. The left inset shows 
formation volume of vacancy relative to the volume of the 
unit cell, while the inset on the right shows the experimental 
supersolid fraction from Ref. Q- 



wiat)- Not only the peak location, but also the overall 
shape of tt vac (P) strikingly resembles the pressure de- 
pendence of the supersolid fraction as measured by Kim 
and Chan [5[ (shown for comparison on the right inset of 
Fig. S]). While we do not have a rigorous explanation for 
this coincidence, it is worth noting that the peak value 
of $7 vac is close to t>i a t. Consider some volume V = L 3 
surrounding a vacancy. To ensure the correct value of 
^vac — flat i the deformation of lattice lines on the border 
of this imaginary volume has to scale as (n vac — vi a t)/L 2 . 
Therefore, the strain on the border of this volume is pro- 
portional to (f2 vac — v\ at )vllt /V. Volume encompassing 
a region with the strain above some threshold e c is given 
by V c ex (£l vaiC — v\ at )v^ / e c . This means that the volume 
of the strained region of the lattice that accompanies a 
vacant lattice site in solid helium is strongly reduced at 
pressures close to 60 bar. 

To conclude, we have been able to characterize the iso- 
baric vacancy formation in solid 4 He. Thermodynamic 
vacancy properties were extracted from the calculations 
of the equations of state of solid helium with number of 
atoms both commensurate and incommensurate with the 
number of available atomic sites. Isobaric vacancy for- 
mation energy turned out to be significantly lower than 
the energy of formation at fixed density. The value of the 
isobaric formation energy is nonetheless high enough to 
exclude the possibility of intrinsic non-interacting vacan- 
cies in the supersolid experiments. However, we find a 
strong resemblance between the pressure dependence of 
vacancy formation volume and experimental results for 
the supersolid fraction. The maximum in the formation 
volume is accompanied by a formation energy minimum 
at the same pressure. These coincidences suggest that 



vacancies may be in fact in some way involved in the 
supersolidity mechanisms. 
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